# 此程序可将WGS84坐标系下的平面坐标值转换为经纬度
# 用于验证录制的坐标值是否真的为WGS84坐标值


from pyproj import Proj, transform
 
def utm_to_latlon(easting, northing, zone_number, zone_letter):
    # 构建UTM投影的Proj对象
    utm_proj = Proj(f"+proj=utm +zone={zone_number} +datum=WGS84 +units=m +no_defs +south")
    if zone_letter.lower() in ['s', 'south']:
        utm_proj = Proj(f"+proj=utm +zone={zone_number} +datum=WGS84 +units=m +no_defs +south")
    else:
        utm_proj = Proj(f"+proj=utm +zone={zone_number} +datum=WGS84 +units=m +no_defs")
    
    # 构建WGS84的Proj对象
    wgs84_proj = Proj(proj='latlong', datum='WGS84')
    
    # 使用transform函数进行坐标转换
    lon, lat = transform(utm_proj, wgs84_proj, easting, northing)
    return lat, lon
 
# 示例：将UTM坐标转换为经纬度
# easting = 703337.858  # 东向距离（米） x 实例数据 南京地图录制原始数据第一行 
# northing = 3482675.900  # 北向距离（米）y  3998395.095,773118.408

easting = 773118.408  # 东向距离（米） x 实例数据 安阳地图录制原始数据第一行 
northing = 3998395.095  # 北向距离（米）y 

zone_number = 49  # UTM区域号  带号=（经度整数位 / 6）的整数部分 + 31  这里要按照手持惯导手簿配置的实际UTM带号进行设置
zone_letter = 'N'  # UTM区域字母（N表示北半球，S表示南半球）
 
lat, lon = utm_to_latlon(easting, northing, zone_number, zone_letter)
print(f"Latitude: {lat}, Longitude: {lon}")